clc; clearvars -except NumberConsecutiveshocks FromM2Q LogVsLevel
  LogVsLevel = 0; 
  frequency = 12; % Monthly
% frequency =  4; % Quarterly
%-------------------------------------------------------------------------%
[x_, ~] =xlsread('./US_Historical_EPU_data.xlsx');      
ss_Bloom = find(x_(:,1)==1960);                                   %se_Bloom = find(x_(:,1)==1984); 
dates_EPU_HIST    = x_(ss_Bloom(12/frequency):12/frequency:end,1); se_Bloom = find(dates_EPU_HIST==1984); 
EconomicPolicyUnc = x_(ss_Bloom(12/frequency):12/frequency:end,3); clear x_ % 1960/01-2014/10
  
[x_, ~] =xlsread('./US_Policy_Uncertainty_Data.xlsx');  
  EconomicPolicyUnc_bis     = x_( 12/frequency:12/frequency:end,3);%clear x_  % Baseline_overall_index
  EconomicPolicyUnc_NewOnly = x_( 12/frequency:12/frequency:end,4); clear x_  % News_Based_Policy_Uncert_Index - to interpret the unit of increase (R2, 2nd round)
  
EPU           = [ EconomicPolicyUnc(1:se_Bloom(end));                                                                             EconomicPolicyUnc_bis];
EPU_NewsOnly  = [ EconomicPolicyUnc(1:se_Bloom(end));                                                                             EconomicPolicyUnc_NewOnly];
EPU_rescaled  = [(EconomicPolicyUnc(1:se_Bloom(end)))*std(EconomicPolicyUnc_bis(1:length(EconomicPolicyUnc(se_Bloom(end)+1:end))))...
                                                                                    /std(EconomicPolicyUnc(se_Bloom(end)+1:end)); EconomicPolicyUnc_bis]; % 1960/01-2019/10

cstrEPU =  cal(1960,1,frequency,length(EPU));
REPLICATION_BBD = 0;
if REPLICATION_BBD
    ssEPU = ical(1985,1,cstrEPU); seEPU = ical(2014,12,cstrEPU);
else
   %ssEPU = ical(1960,1,cstrEPU); seEPU = ical(2019,12,cstrEPU);
    ssEPU = ical(1985,1,cstrEPU); seEPU = ical(2019,12,cstrEPU); % only NEW EPU available since 1985/01
end
if LogVsLevel
  logEPU          = log(EPU(         ssEPU:seEPU));
  logEPUnews      = log(EPU_NewsOnly(ssEPU:seEPU));
  logEPU_rescaled = log(EPU_rescaled(ssEPU:seEPU));
else
  logEPU          =     EPU(         ssEPU:seEPU) ; 
  logEPUnews      =     EPU_NewsOnly(ssEPU:seEPU) ; 
  logEPU_rescaled =     EPU_rescaled(ssEPU:seEPU) ;
end

clearvars -except NumberConsecutiveshocks logEPU logEPU_rescaled logEPUnews REPLICATION_BBD FromM2Q

%% From VAR as in BBD QJE
[x_, ~] =xlsread('./BBD_VAR_DATA.xlsx');   
if REPLICATION_BBD
    ss_ = find(x_(:,1)==198501);  se_ = find(x_(:,1)==201412); 
else
   %ss_ = find(x_(:,1)==196001);  se_ = find(x_(:,1)==201912);
    ss_ = find(x_(:,1)==198501);  se_ = find(x_(:,1)==201912);
end

%dataMatrix= 100*[logEPUnews/100       log(x_(ss_:se_,2)) x_(ss_:se_,3) log(x_(ss_:se_,4:5))]; 
 dataMatrix= 100*[logEPU_rescaled/100  log(x_(ss_:se_,2)) x_(ss_:se_,3) log(x_(ss_:se_,4:5))]; 
dataName   = {'EPU','SP500 index','Short-term nominal rate','Employment', 'Industrial production'}; 
clear ss_ se_ x_
VarsToPlot = [5 4]; pos_shock = find(strcmp(dataName,'EPU'));% pos_shock = 1;

nobs       = size(dataMatrix,1);

%-------------------------------------------------------------
% Step 0: parameter settings
%-------------------------------------------------------------
nlags  = 6;                   % number of lags in VAR  
nconst = 1;                   % 3 = for constant, linear and quadratic trend, 2 = for constant and trend, 1 = for constant, 0 otherwise 

%-------------------------------------------------------------------------%
prior=1;                      % for prior,  1   = standard prior
                              %             inf =     flat prior

%set prior for first lag of each variable (if estimation with Minnesota prior)
[~,coly]=size(dataMatrix);
lev=ones(coly,1);                                          
%-------------------------------------------------------------------------%       
  nimp = 36;                  % horizon of VD and IRFs to be shown         

est_method = 0;               % 0 : OLS point estimation
                              % 1 : Bayesian estimation with Minnesota prior 
                              % 2 : OLS point estimation and classical bootstrap    
                              
%-----------------------------------------------------------------------
% Step 1: load and transform data
%-----------------------------------------------------------------------
[y,x] = mk_vdm(dataMatrix,nlags,0);
%y is T x nvars matrix of rhs variables
%x is T x nvars*nlags  of lhs variables
[T,nvars]=size(y);

%---------------------------------------------------------------------
% Step 3: Estimate VAR, extract shocks, FEV shares explained, and IRFs 
%---------------------------------------------------------------------
if nconst==1
   const = ones(T,1); 
   x = [x const];
   ncoeffs = nvars*nlags+1;
elseif nconst==2
   x = [x ones(T,1) (1:T)'];
   ncoeffs = nvars*nlags+2;
elseif nconst==3
    x = [x ones(T,1) (1:T)' ((1:T).^2)'];
%   dummy_1975Q2 = dates_q_==1975.5;
%   x = [x ones(T,1) (1:T)' ((1:T).^2)' dummy_1975Q2(1:end-nlags)];
   ncoeffs = nvars*nlags+3;   
else
   ncoeffs = nvars*nlags; 
end

%ncoeffs + 1 x nvars matrix of regression coefficients
b=x\y;       %OLS coefficients: ncoeffs x nvars (** more efficient way of computing least squares; i.e. b=inv(x'*x)*x'*y;
res = y - x*b;          %T x nvar matrix of residuals
vmat = (1/T)*(res'*res);
stds=sqrt(diag(vmat));

%OLS point estimates
[FEVD_IRF,v]  =RecursivenessIdentification(b,res,vmat,nvars,nlags,nimp, pos_shock);

%----------------------------------------------------------------------
% Step 3: Plots
%----------------------------------------------------------------------
% clearvars -except dataName VarsToPlot FEVD_IRF plottype nvars
% save(strcat('IRFDATA_',plottype))

% %FEV shares explained by shock
% plotFEV(FEVD_IRF,vars,slope,use_yields,shock_extract);
% 
%IRFs to shock
% plotIRF(FEVD_IRF,vars,slope,use_yields,shock_extract)

if REPLICATION_BBD
% FEVD_IRF
[nimp,~,~]=size(FEVD_IRF);
% R=3; %round(N/2);
h=figure('Color',[0.9412 0.9412 0.9412],'Position',[1 1 800-100 600-100],'Name','IRFs to 1st shock');
figure(h);
plotnum=0;
for n=VarsToPlot
    plotnum=1+plotnum;
    if     length(VarsToPlot)==2 
        subplot(1,2,plotnum);             
    elseif length(VarsToPlot)==5
        subplot(2,3,plotnum);   
    end

    if ndims(FEVD_IRF)==3   %if available, plot confidence bands
        grpyat = [(1:nimp)' FEVD_IRF(1:nimp,nvars+n,2); (nimp:-1:1)' FEVD_IRF(nimp:-1:1,nvars+n,3)];
        patch(grpyat(:,1),grpyat(:,2),[0.7 0.7 0.7],'edgecolor',[0.65 0.65 0.65]);
    end
%     end
    hold on;
% /// responses of industrial production and employment to a 90-point upward EPU innovation,
% ~ 4 std dev
%   plot(1:nimp,       FEVD_IRF(:,nvars+n,1)                      ,'k-');%,'LineWidth',2);
    plot(1:nimp,3.7292*FEVD_IRF(:,nvars+n,1)                      ,'k-');%,'LineWidth',2);
%     end
%   plot(1:nimp,FEVD_IRF(:,nvars+n,1)/FEVD_IRF(1,nvars+2,1),'k-');%,'LineWidth',2);
    plot(1:nimp,zeros(nimp),':k');
    title(dataName(n),'FontSize',14) %'FontWeight','bold',
    ylabel('percent' ,'FontSize',12)
    xlabel('Months','FontSize',12)
    axis tight
   %axis([0 nimp ylow(n) yhigh(n) ]);
    hold off;
end
end
uEPUVAR = v(12-nlags+1:end);
clearvars -except NumberConsecutiveshocks logEPU logEPUnews logEPU_rescaled uEPUVAR FromM2Q

%% FIT AR(p) as in Brogaard and Detzel (MS, 2015)
sysOP = 'MAC';
if strcmp(sysOP,'MAC')
    x_ =  readtable('./F-F_Research_Data_Factors.CSV');      
else
   [x_, ~] =xlsread('./F-F_Research_Data_Factors.CSV'); 
end
%ss_ = find(x_{:,1}==196001);
 ss_ = find(x_{:,1}==198501);
 se_ = find(x_{:,1}==201912); 
dates_EPU = x_{ss_:se_,1};clear x_ ss_ se_% 1960/01-2014/10
%--> Shocks 1961/M1 -- 2019/M12
dates_EPU = dates_EPU(13:end);
%%
    uEPU = uEPUVAR  ; clear  uEPUAR1 uEPUAR12           uEPUVAR 
    
    uEPU = uEPU/std(uEPU);
% /////////////////////////////////////////////////////////////////////////
%----------- Construct indicators for seq. of shocks with 2 ident schemes-%
IndPreviousNShocksEPU= zeros(size(uEPU));
% /////////////////////////////////////////////////////////////////////////
for tt=NumberConsecutiveshocks:length(uEPU)-1 
    if tt==NumberConsecutiveshocks
%      %if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPU(tt+1)<0);  IndPreviousNShocksEPU(tt) = 1; end
        if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksEPU(tt) = 1; end
    else
%      %if (uEPU(tt-NumberConsecutiveshocks)<0 && sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPU(tt+1)<0);  IndPreviousNShocksEPU(tt) = 1; end
%       if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksEPU(tt) = 1; end
        if (uEPU(tt-NumberConsecutiveshocks)<0 && sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksEPU(tt) = 1; end
    end
end
clear tt

logEPU          = logEPU(         13:end);
logEPU_rescaled = logEPU_rescaled(13:end);